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Abstract The effect of structural disorder on the stress re- 
sponse inside three dimensional particle assemblies is stud- 
ied using computer simulations of frictionless sphere pack- 
ings. Upon applying a localised, perturbative force within 
the packings, the resulting Green 's function response is mapped 
inside the different assemblies, thus providing an explicit 
view as to how the imposed perturbation is transmitted through 
the packing. In weakly disordered arrays, the resulting trans- 
mission of forces is of the double-peak variety, but with peak 
widths scaling linearly with distance from the source of the 
perturbation. This behaviour is consistent with an anisotropic 
elasticity response profile. Increasing the disorder distorts 
the response function until a single-peak response is ob- 
tained for fully disordered packings consistent with an isotropic 
description. 

Keywords Force Perturbations • Stress Response 
PACS 61.43.-j • 81.70.Bt • 83.80.Fg 



1 Introduction 

The internal stresses of a particulate medium, such as a sand- 
pile, determine mechanical stability and structural robust- 
ness. How these stresses are distributed can have a profound 
effect on the way such a system responds to external load- 
ing. This is a challenging problem which, one hopes, will 
lead to a better understanding of the transition between qui- 
escence and catastrophic failure, such as a collapsing grain 
silo or the onset of an avalanche. 

Over the past decade or so, it has come to light that 
the way forces are transmitted through a disordered pile of 
grains occurs in a rather inhomogeneous manner B23I1 . It 
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seems that it is very much the particulate nature of a gran- 
ular material that lends itself to the emergence of strongly 
heterogeneous force networks in disordered packings. Con- 
sequently the nature of stress transmission in response to 
external loading might also exhibit unusual properties due 
to granularity. These issues have led to renewed interest in 
describing how stresses are transmitted through particulate 
matter. 

Whereas classical, isotropic, elasticity theory EEl 
accurately describes the stress state of continuous media, 
when the discrete nature of the constituent particles influ- 
ence the way a system responds, it is unclear how one should 
proceed. So much so that recent simulations of Lennard- 
Jones glasses have shown how the discrete nature of the 
amorphous glass influences the elastic properties of the ma- 
terial Jill [la I27I Est - Thus, in the context of jammed partic- 
ulate media, the humble sandpile has thus come to represent 
the paradigm of an unusual state of matter exhibiting hetero- 
geneous force distributions. 

One of the simplest methods to access information about 
the stress state of a material is to measure the Green's func- 
tion response to a localised perturbation. For a granular pack- 
ing, this would involve applying a small force to the central 
particle at the top of the packing then measuring the pressure 
at the bottom of the packing in response to this localised ex- 
ternal perturbation. To remain in the linear regime, the ap- 
plied force perturbation must be sufficiently small so as to 
avoid rearrangement of the grains from their original posi- 
tions. This is the approach taken in several experimental and 
simulation studies to date. 

Some of the more striking studies involved the direct vi- 
sualisation of the force networks in two dimensional (2D) 
packings using photoelastic particles HUSHED, an d their 
micro-displacements |0;E3l- Direct imaging and the sub- 
sequent analyses allows one to determine how the contact 
and/or force network responds when a point force is ap- 
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plied to the packing. Presently, more indirect methods are 
required for inferring the way forces are transmitted through 
3D packings Q Qfl H H, simply due to the fact that 
granular packs are inconveniently opaque. However, there is 
a clear distinction between the response properties of dis- 
ordered grain piles [33] and ordered arrays [30]. There are 
also numerous analytical and numerical works 111 ; 12 ; 13 
SlilllliiliilllllllSllillliilillii, predicting how 



forces are transmitted in response to such perturbations. 

As a precis to the above, also see Ref. [20], it is now ap- 
preciated that granular materials have the uncanny ability to 
respond in a variety of ways. Isotropic, single-peak response 
profiles mimic the manner in which an elastic medium be- 
haves as described by partial differential equations of the 
elliptic class. The isotropic elastic stress state, a, of a ma- 
terial in response to a localised load, F, on an infinite half 
space is described by the Boussinesq equation [24], 
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where z is the vertical distance, or depth, from the source of 
the force perturbation, and r is the planar distance from the 
axis of F 1126c 12411 . Most of the experiments and simulations 
performed to date are constructed in such a manner as to di- 
rectly test the validity of this theory and its 2D analogue the 
Flamant equation. Under some conditions, strongly anisotropic, 
double-peakresponse profiles have been reported. These types 
of response profiles can be indicative of anisotropic, elas- 
tic behaviour, and therefore formally belong to an elliptic 
description. In the extreme case, double-peak behaviour is 
representative of hyperbolic, wave-like propagative models, 
whereby the stresses propagate along characteristic path- 
ways analogous to light rays. In between, is a diffusive parabolic 
model that seems to be losing favour despite playing a sem- 
inal role in describing force heterogeneities in granular me- 
dia. Moreover, a crossover between different response func- 
tions can occur, even within the same pile, depending on 
structure, particle friction coefficient, distance from the source 
of the perturbation, and the magnitude of the imposed force. 
Typically, however, ordered grain configurations result in 
anisotropic response features, whereas disordered systems 
usually conform to an isotropic picture. 

As far as 2D results go, many of the features described 
above can be accommodated within the framework of 2D 
anisotropic elasticity theory lfl7[ l32ll . In fact, a parameter 
space of anisotropic elasticity theory was proposed IB2I1 that 
allows for the appearance of either single-peak or double- 
peak response functions within the same formalism. Sim- 
ply put, anisotropy in the force/contact network promotes 
an anisotropic response. One of the key concerns, therefore, 
is the ability to distinguish between the elastic and wave-like 
response profiles. 



Clearly there are many factors that influence the nature 
of the resulting response. It would seem prudent to investi- 
gate the underlying factors that have the greatest influence 
in determining the response properties. This is precisely the 
approach taken in this study. Structure, or the arrangement 
of the particles that constitute the packing, plays a dominant 
role. The results presented here are the first to investigate the 
manner in which the stresses inside 3D particle packings are 
transmitted in response to localised force perturbations, and 
how varying the disorder affects the response. 



2 Simulation Model 

The computer experiments reported here are designed to sys- 
tematically study how structure influences stress transmis- 
sion in response to localised perturbations. To that end, the 
simulations are carried out on a model system: three dimen- 
sional, monodisperse, non-cohesive, frictionless, soft-sphere 
packings, with fully periodic boundary conditions, in the 
absence of gravity. Generation of the initial packings has 
been described in detail elsewhere 1134c 13511 . To summarise, 
N spheres of diameter d and mass m, were arranged into a 
face-centred cubic (fee) array. The packing fraction, 0, of 
the particle assemblies was fixed at (j) — 0.742, just above 
that of a hard-sphere fee array, 0f cc = ^/2k/6. Two parti- 
cles, i and j, are defined to be contact neighbours and in- 
teract through a short-range, purely repulsive, force when 
ry < d, where ry is their centre-centre separation. The par- 
ticles interact via a Hookean force law, = k(d — ry), for 



ry < d, and zero otherwise, i.e. a one-sided linear spring. In 
this study, the particle stiffness k merely parameterises the 
force scale and for convenience is set to unity, along with 
d and m. Distances are reported in units of d and all other 
quantities are appropriately non-dimensionalised. 

Configurations with different amounts of disorder were 
generated by adding defects to the original fee array. In this 
study, these defects were introduced by randomly removing 
particles from the original array then allowing the configu- 
ration to relax into a different mechanically stable state at 
the same 0. Removing more particles prior to the relaxation 
process resulted in more disorder. As a result the number of 
particles varied from N = 16384 for the FCC lattice down 
to N = 10384 for the most random packings studied. Con- 
sequently, the size of the simulation cube, L, ranged from 
19.4c/ < L < 22.6c/ The amount of disorder was quantified 
by the coordination number z (average number of contacting 
neighbours), given in Table [U which varied from 12 for the 
fee array, to below 9 for the most disordered packing. Other 
measures of the packings corroborated the gradual transi- 
tion from ordered to disordered configurations j34cl35ll . This 
configuration-generation protocol provided a suitable man- 
ner in which to control the amount of disorder. The pack- 
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Table 1 Configurations are labelled by Ci and their coordination num- 
ber z. The fee array has a value = 12. 



Ci 


FCC 


C2 


C4 


C6 


C8 


CIO 


C12 


z 


12 


11.85 


11.71 


11.22 


10.68 


10.29 


8.92 



ings are labelled accordingly: C[\ —6] correspond to quasi- 
ordered assemblies, C[7— 11] weakly disordered, and C12 
fully disordered [35]. Table Q] catalogues the packings used 
in this work with their respective coordinations. These data 
were obtained after averaging over five independent realisa- 
tions. 

The major aim of this work is to investigate how pack- 
ing structure affects the way localised force perturbations 
are transmitted through sphere packings. In an effort to elim- 
inate other contributing factors that may strongly influence 
the response properties, a particular geometry is used to study 
stress transmission. Specifically, to probe the nature of stress 
transmission inside the packings, a localised, isotropic per- 
turbation was imposed within the central region of the pack- 
ing by increasing the diameter of the centrally located par- 
ticles by an amount Sd, in the range 10~ 4 < Sd/d < 10~ 2 . 
These values for 8d are directly converted into forces using, 
Fpenurb = k(Sd/2). Although this method of perturbation is 
not standard compared to experimental and previous 2D nu- 
merical studies, it provides a convenient method to study the 
response properties of unconfined packings in the absence 
of sidewalls and directional external forces such as gravity. 
Moreover, this method also improves the statistical analysis 
as discussed below. 

For the results presented here to connect with expected 
theories and experimental works, i.e. to satisfy linearity, it is 
necessary to gauge the range of the magnitude of the pertur- 
bative force that can be applied to the packing. This value 
was based on the average force per particle, f av , in the static 
packings prior to the perturbation procedure. These average 
forces varied over, 10~ 3 ta/ < / av < 10~ l kd, from the or- 
dered to most disordered configurations. Thus, the strength 
of the perturbations was varied from below to above the av- 
erage force in most cases. After the application of this per- 
turbation, the system was again relaxed at the same initial 
= 0.742. This allowed a direct comparison between the 
initial stress state a 1 of the configuration before the pertur- 
bation and the stress state <7 f , of the final configuration af- 
ter. The contact stresses a were computed for all contacting 
particles inside a volume £2 (l]], o a p = j2T.ij r fjf!j^ where 
the a,j3 are the cartesian components of particles i and /', 
separated by ry. The effect of the imposed perturbation was 
computed as the change in the normal stress between the 
final and initial states, a = <J f — &. To reduce statistical 
uncertainty, for each independent realisation, the stress was 
averaged over the equivalent faces of the periodic cell, i.e. 
° = + °» + °vy + °yy + °u. + °a)A>» where the + and 



— represent the normal stresses calculated in the positive and 
negative direction relative to the position of the perturbation 
source. See Fig.Q] 




Fig. 1 Schematic of the simulation configuration. The perturbed region 
is identified by the particle in the middle of the packing. The planes of 
regular sized spheres represent examples of the visualisation planes 
for the stress transmission presented in the following figures. All other 
spheres are drawn at a reduced size for clarity. Axes are labelled al- 
though they are all equivalent in this geometry. The simulation cube 
has side length L = 22.6d. 

Figure Q] presents a schematic of the simulation geom- 
etry. The shaded regions represent the planes over which 
the stresses were computed and visualised in the following 
figures. The dark gray plane is denoted the direction paral- 
lel to the direction of stress transmission, whilst the light 
gray plane is perpendicular to the transmission direction. 
The stress response in the plane above the (+z) and below 
(— z) the point of perturbation are equivalent because, a) the 
packings are isotropically compressed, b) there is no gravity 
acting on the system (i.e. directional external forces), and c) 
because the perturbation itself is isotropic. Likewise, the +x 
and — x and +y and — y are also all equivalent. Observing the 
stresses in these planes allows a direct comparison to 2D re- 
sults. This procedure might appear somewhat artificial, how- 
ever, there is a close resemblance between this simulation 
protocol and previous experimental studies. One must sim- 
ply imagine taking the simulation cube and placing any of its 
faces on a solid floor and then measuring the pressure pro- 
file at this boundary. The benefit of this particular simulation 
scheme is the ease with which one can gain 'depth' depen- 
dent information, where 'depth' is the mean distance from 
the source of the localised perturbation. Therefore, rather 
than generating different packings at different heights, one 
can simply look at different slices within the packing J2]. 
Moreover, the protocol implemented here could be realised 
in experiments on colloidal glasses and foams. 
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3 Results 

Figure [2] shows stress response maps for the weakly disor- 
dered packing C2, for one value of the perturbation force. 
The top panels correspond to different planes perpendicular 
to the direction of stress transmission (c.f. light gray plane in 
Fig-ID- The different panels of Fig.|2]show the stress at dif- 
ferent 'depths' h, with h = defining the central plane of the 
packing, coincident with the perturbation source. While the 
response function gradually decays with distance from the 
perturbation, in the immediate vicinity of the perturbation 
(top left panel), the response is localised (single lobe), cap- 
turing the finite extent of the imposed perturbation. This dra- 
matically changes character away from the source of the per- 
turbation where the response becomes strongly anisotropic 
(ringed response). The bottom panel in Fig. [2] shows the 
plane parallel to the direction of stress transmission (c.f. 
dark gray plane in Fig. [TJ for the same C2 packing. The 
response function is directed along characteristic, broaden- 
ing pathways, leaving the region between these high-stress 
paths experiencing a relative depletion of stress. This view 
is similar to the picture put forward in Ref. H30I1 on stress 
transmission through ordered arrays of glass beads. 




Fig. 2 Stress response inside a weakly disordered packing, C2, due to 
perturbative force Sd = 0.0001. Top panel: Stress maps of approxi- 
mately ~ 22.6d x 22. 6d in size, corresponding to three planes parallel 
to the light gray plane in Fig.fTJ respectively located at h = 2\/2, 3\/2, 
and 5 \f2, from left to right. There is a qualitative change in the nature 
of the response with increasing distance from the perturbation source. 
The strength of the response also decays with distance from the source. 
Bottom: Stress map for a central plane corresponding to the dark gray 
plane in Fig.fT] with dimensions ~ 22.6d wide along the x axis and ex- 
tending out ~ I2d in the z direction away from the perturbation source. 
Darker shades indicate larger a. Intensity scale same in all panels. 



To further quantify the nature of the response, and in 
an effort to compare with equivalent 2D results, the stress 



profiles shown in Fig. [3] are central cuts taken from {e.g., 
the top panels in Fig. |2]i slices at different h (equivalently, 
e.g. bottom panel in Fig. |2j inside the C2 packing for differ- 
ent perturbation magnitudes 8d. This clearly demonstrates 
that close to the source the response is single-peaked, cross- 
ing over to a double-peaked response further away. Not only 
are the rescaled stresses independent of 8d, but no contacts 
were broken in the perturbation process, thus confirming the 
linear response regime. To check the robustness of these re- 
sults, qualitatively similar data were obtained when k was 
varied over 5 orders of magnitude, and the degree of over- 
compression, (j) — 0f cc , covered 3 orders of magnitude (re- 
sults not shown here). 




Fig. 3 Stress profiles rescaled by the 'depth', h, from perturbation 
source and the magnitude of the force, Sd, taken from central cuts 
across the planes shown in the top panels of Fig.f2](packing C2). The x- 
axis denotes the distance across the packing face. Each trio of data are 
for h = 2\/2 (circles), 3\/2 (squares), and 5\/2 (diamonds), for differ- 
ent values of the perturbing force: Sd = 0.0001 (open), 0.001 (filled), 
and 0.01 (grey line). 



The results of Figs.|2]and[3]indicate that away from the 
source of the perturbation, the response functions point to- 
ward possible agreement with the hyperbolic wave theory 
of force propagation II 1011 . which predicts a double-peaked 
response function with peak widths, W, scaling diffusively 
with depth, W ~ \[h. To determine which class of theory 
these results belong to, Fig. [4] shows the depth dependence 
of W and the peak separations 4r pea k. From Fig.|4]the char- 
acteristic features of the double-peaked response function 
are seen to vary linearly with depth, W ~ h, and therefore 
do not belong to the hyperbolic class of solutions. This sug- 
gests that the response properties of the weakly disordered, 
frictionless arrays studied here are likely to be described by 
three dimensional anisotropic elasticity theory. 

Although a developed formalism for anisotropic elas- 
ticity exists, there are no analytic solutions with which to 
directly compare these simulation results. Alternatively, an 
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line) is shown. A criterion was used in the fitting procedure 
that attempted to match the peak positions as closely as pos- 
sible between the simulation data and Eq|2] Despite these 
minor inadequacies, Eq. |2] can be used to characterise the 
profiles in the weakly disordered regime through the coef- 
ficients c' and co, which are obtained as fitting parameters. 
Over the range of disorder explored in Fig. Ha), the values 
of these fitting parameters are shown in Fig. Qb). The up- 
per left points correspond to the strongly ordered systems 
whereas the lower right points correspond to the weakly dis- 
ordered configurations. These coefficients not only follow 
the same trend, but are somewhat similar, to those obtained 
in ordered and weakly-disordered 2D experiments llal . 



Fig. 4 Stress peak analysis for configuration C2 that exhibits a double- 
peak response function. The width W (•), measured at half height of the 
stress response peak averaged over the two peaks, and the separation 
between the two peaks, A (o). 



empirical fitting form can be used to describe the data. Moti- 
vated by twin-peak, hyperbolic models, Breton et al. |@] in- 
troduced the Convection-Diffusion (CD) equation that mim- 
ics the response properties of a 2D lattice. The CD equation 
is characterised by a propagation, or 'speed', coefficient c, 
that determines the direction of the stress response, and a 
'diffusion coefficient' D, characterising the diffusive spread- 
ing of the two peaks with depth. In its original form, the 
resulting scalings have the width of the peaks, W — \j2Dh, 
consistent with a wave-like hyperbolic picture. This was mod- 
ified into the Convection-Wave (CW) equation by Geng et al. 
1 1511 . in an effort to deduce whether their experimental data 
was better described by the hyperbolic approach or more in 
line with an elliptic model. To that end, the CW equation is 
extended here for 3D systems, 

F 



a(h) = 



2(2nff 2 aP-h 



-Al/lm 2 !? , -A\l2m 2 h 2 



(2) 



F is the magnitude of the perturbing force, r is the in-plane 
radial distance from the central axis of the packing, and 
A± = (r ± c'h). Equation|2]guarantees a twin peak response 
with peak widths scaling linearly with depth, W = coh, where 
CO and c' depend on the packing structure. 

Figure[5Ja) shows the stress response profiles at one depth 
for a given perturbation, over a range of weak disorder. The 
double-peak profiles evolve with disorder: the amplitude of 
the response decreases with increasing disorder, although 
over this range in disorder, there appears to be little change 
in the location of the peaks. Equation [2] captures the main 
features of the profiles - the characteristic double-peak - yet 
it fails to accurately describe the full response, such as the 
dip in the profiles either side of the peaks. A representative 
fit of Eq. |2] (symbols) to the C6 configuration (thick solid 
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Fig. 5 (a) Scaled stress profiles at h = 3\/2 for 5 = 0.01, in the 
weakly disordered regime. Peak heights decrease with increasing dis- 
order from highest, C2, to lowest, C8. The symbols (+) are a fit of Eq.[2] 
to the data for C6 (thick solid line), with best-fit parameters for c' and 
CO, and is representative of the goodness of the fit for the other configu- 
rations, (b) The values of the fitting parameters c' and CO used in Eq.ff] 
extracted from fitting to the packings C2 — C8. 
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When further disorder is introduced, the anisotropic re- 
sponse profiles become less well-defined and there is a crosso 
to a more isotropic response of the single-peak variety. This 
is illustrated in the stress maps shown in Fig.|6l for configu- 
rations with moderate-large disorder. In an effort to improve 
averaging over the growing stress fluctuations, both the size 
of the perturbed region and the averaging volume were var- 
ied. It is expected that even for frictionless particles, disor- 
der promotes an isotropic, elastic-like response 11911 . For the 
most disordered system (CI 2), whose structure resembles 
that of an overcompressed random closed packed state, a 
single-peaked response can be clearly resolved. This is bet- 
ter illustrated in the profile shown in Fig. [7] A comparison to 
the Boussinesq theory of Eq.Q] shows reasonable agreement 
fl. Further promoting the idea that isotropic random pack- 
ings are more likely to exhibit isotropic elastic-like stress 
behaviour. 




Fig. 6 Stress response at h = 3v2 for 8d = 0.001, inside packings 
with increasing disorder from left to right: C8, CI 1, and C12. Size of 
slice ~ 2Qd x 20d. Shading is the same in all panels. 



4 Discussion 

The procedure described here presents only one of a num- 
ber of avenues through which the mechanical response of a 
material may be investigated. Traditionally, shear tests have 
been used to determine constitutive parameters through which 
material, or elastic, constants can be obtained from linear 
stress-strain relations. Moreover, for the explicit study of re- 
sponse properties of the type studied here, the elastic con- 
stants must be initially determined which can then be input 
into finite element analyses routines |4,]. However, in terms 
of computational time such procedures are just as costly, if 
not more so, as the perturbation protocol implemented here. 

To reiterate, the goal here was to make a connection 
with experimental studies on similar systems. The results 
presented here are the first in a series of investigations on 
the mechanical response in particulate media. In this initial 
study it has been shown that the force response method pro- 
vides a quantitative description of the elastic properties of 
particle packings consistent with previous studies. To that 
end, the protocol implemented here has been verified as a 
suitable procedure through which we can investigate me- 
chanical response problems. Even though the particular ge- 
ometry used is not conventional from an experimental stand- 
point it does offer a method to determine the influence of 
different parameters on the mechanical behaviour of sphere 
packings. The focus here was structure and is likely rel- 
evant not only to compressed elastic sphere-packings, but 
also foam dispersions. 
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Fig. 7 Comparison between the Boussinesq theory, Eq. [T] (dashed 
line) and the stress response for the most disordered configuration 
CI 2, which resembles an overcompressed random close packing. At 
h = 3 </2 inside the packing for S d = 0.001. 



5 Conclusions 

For the first time, a systematic study of the nature of stress 
transmission in response to localised force perturbations has 
been obtained inside three dimensional packings of friction- 
less spheres using discrete element simulations, in a particu- 
lar geometry that avoids the influence of walls. The manner 
in which stresses are transmitted crucially depends on the 
distance from the source of the perturbation and the struc- 
ture of the packing. For small deviations from the fee struc- 
ture, the response in the vicinity of the perturbation source 
retains the localised form of the perturbation, and appears 
as a single-peak response function. This crosses over to an 
anisotropic response function of the double-peak variety fur- 
ther from the source. This double-peak feature is suitably 
captured by the modified Convection- Wave equation; peak 
widths scaling linearly with distance, rather than the diffu- 
sive form of this equation [9]. These results are consistent 
with a generalisation from 2D to 3D anisotropic elasticity 
theories 03211 . Further disorder deforms the double-peak fea- 
tures associated with the anisotropic response function until 
there is a crossover to a single-peak response for randomly 
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packed spheres consistent with an isotropic elastic descrip- 
tion. Moreover, this crossover is qualitatively matched by 
the Boussinesq equation, Eq. Q] which again indicates that 
the procedure implemented here probes mechanical proper- 
ties that is simple and quite accurate [41]. This will be par- 
ticularly useful in future studies on random packings with 
different coordination numbers with and without friction. 

It is known experimentally the stress response of a disor- 
dered, frictional granular packing qualitatively agrees with 
isotropic elasticity theory ll33ll . Simulations of two dimen- 
sional systems also point towards agreement with isotropic 
elastic models [19]. Studies of three dimensional frictional 
ordered and disordered packings and gravity packings re- 
mains work in progress. 
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